Penalty Robin-Robin domain decomposition 
schemes for contact problems of nonlinear elastic 

bodies 

Ihor I. Prokopyshyn * Ivan I. Dyyak ' 
Rostyslav M. Martynyak * Ivan A. Prokopyshyn § 

September 7, 2012 



Abstract 

In this paper we propose on continuous level several domain decompo- 
sition methods to solve unilateral and ideal multibody contact problems 
of nonlinear elasticity. We also present theorems about convergence of 
these methods. 
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1 Introduction 

Many domain decomposition techniques for contact problems have been pro- 
posed on discrete level, particularly substructuring and FETI methods IlL m ■ 
Domain decomposition methods (DDMs), presented in 0, E3, U, Ell for 



unilateral two-body contact problems of linear elasticity, are obtained on con- 
tinuous level. All of them require the solution of nonlinear one-sided contact 
problems for one or both of the bodies in each iteration. 

In works [II [01 El we have proposed a class of penalty parallel Robin-Robin 



domain decomposition schemes for unilateral multibody contact problems of 
linear elasticity, which are based on penalty method and iterative methods for 
nonlinear variational equations. In each iteration of these schemes we have to 
solve in a parallel way some linear variational equations in subdomains. 

In this contribution we generalize domain decomposition schemes, proposed 
14, 15, (j[ to the solution of unilateral and ideal contact problems of nonlinear 



elastic bodies. We also present theorems about the convergence of these schemes. 
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2 Formulation of multibody contact problem 



Consider a contact problem of N nonlinear elastic bodies f2 Q C M 3 with sec- 
tionally smooth boundaries T a , a = 1,2, N (Fig. 1). Denote f2 = U«=i ^a- 




Fig. 1 Contact of several bodies 

A stress-strain state in point x = (xi, x 2 , x^) T of each body f2 a is defined by 
the displacement vector u a — u a j ej , the tensor of strains e a — e a ij ej and 
the tensor of stresses a a = a a ij &j . These quantities satisfy Cauchy relations, 
equilibrium equations and nonlinear stress-strain law Q : 

O a ij — ^ a $ij ~t~ 2 fl a S a ij 2 /i Q 0J a (e Q ) B a ij , %, j — 1,2,3, (1) 

where Q = £ a n + £ a 22 + £«33 is the volume strain, A Q (x) > 0, > 

are bounded Lame parameters, e a ij — £ a ij — <5ij @a/3 are the components of 
the strain deviation tensor, e a = \J2 g a /3 is the deformation intensity, g a = 

{e all -e a2 2) 2 + (£q22 ~£a33) 2 + {£a33 ~ Sail Y + 6 {e 2 al2 + £ a23 + £ a31 ) 1 and 

is nonlinear differentiable function, which satisfies the following properties: 

Q<u a {z)<d{zLo a {z))/dz<\, d{uj a {z))/dz > 0. (2) 

On the boundary r a let us introduce the local orthonormal basis £ Q , T] a , n a , 
where n a is the outer unit normal to F Q . Then the vectors of displacements and 
stresses on the boundary can be written in the following way: u a — u a ^^ a + 

^ar] a ~\~ ^an ^ ot ) &a — & <\ ' n a &a£ £a ® ar\ V r, ~t~ ^ an 1 a • 

Suppose that the boundary r Q of each body consists of four disjoint parts: 

r Q = r«|jr£Ur£USa, ^ 0, = r«, ?i\}s a ^ 0, where s a = 

{Jp eBa S a /3, and = Ua'ei ^afi'- Surface S a p is the possible unilateral con- 
tact area of body fl a with body fig, and £? Q C {1,2, N} is the set of the 
indices of all bodies in unilateral contact with body il Q . Surface T a pi = Tp> a is 
the ideal contact area between bodies Sl a and Qp>, and I a C {1, 2, N} is the 
set of the indices of all bodies which have ideal contact with il a . 

We assume that the areas S a p c r Q and Sp a C are sufficiently close 
(S a p ra S^), and n„(x) -Wx'), x € S Q/3 , x' = P(x) € S^, where P(x) 



is the projection of x on 5 Q( g [12(. Let <i Q/ 3(x) = ±||x — x'|| 2 be a distance 
between bodies tt a and f2^ before the deformation. The sign of d a p depends 
on a statement of the problem. 

We consider homogenous Dirichlet boundary conditions on the part , and 
Neumann boundary conditions on the part T^: 

u«(x) = o,xer: ; <r Q (x) = p Q (x), xe r*. (3) 
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On the possible contact areas S a p, ft G B a , a = 1,2, ...,N the following 
nonlinear unilateral contact conditions hold: 

cr ml (x) = 073„(x') < 0, ct q? (x) = CT/jf^x') = cr Q?) (x) = ap v (x-) = 0, (4) 

m q „(x) + M/3„(x') < d 0l g(x) , (5) 

KnW + U/3n(x') - d a ^(x) ) cr Q „(x) = , X G S aj a , X = P(x) £ S^ a . (6) 

On ideal contact areas T a pi = Tp> a , /3' G I a , a = 1, 2, iV we consider ideal 
mechanical contact conditions: 

U a (x) = U /3 '(x) , 0" a (x) = (Tpi (x) , XgT^/. (7) 

3 Penalty variational formulation of the 
problem 

For each body fl a consider Sobolev space V a = [i? 1 (0 Q )] 3 and the closed sub- 
space V® — {u a G V a : u a = on T^}. All values of the elements u Q G V^, 
u Q 6 V^ on the parts of boundary T a should be understood as traces [9j . 

Define Hilbert space Vo = V]° x ... x V$ with the scalar product (u,v)y = 

^2^=1 ( u "' v a)y„ and norm ||u||y o = ^J(u , u)y o , u,v G Vo- Introduce the closed 
convex set of all displacements in Voj which satisfy nonpenentration contact 
conditions ([5]) and ideal kinematic contact conditions: 

K = {u G Vo : «„ + ^„<4/i on S a p, u a / = up> on r a ^- } , (8) 

where {a, /?} G Q, Q = {{a, /?} : a G {1, 2, N} , /3 G B Q }, {a', /?'} G Q J , 
Q 1 = {{a', (3'} : a' G {1, 2, iV} , /3' G /„}, and d aP G ff^ 2 (3a), S Q = int (T a \ 

Let us introduce bilinear form A(u, v) = Yl a =i a a(ua,v ), u, v G Vo, which 
represents the total elastic deformation energy of the system of bodies, linear 
form L (v) = ^ a=1 £ Q (v a ), v G Vo, which is equal to the external forces work, 
and nonquadratic functional H (v) = Yl a —i h a (v a ), v G Vo, which represents 
the total nonlinear deformation energy: 

a Q (u Q ,v Q )= / [ X a Q a (u a ) Q (v (V a ) ] t«l, (9) 

■to- « 

ia(v a ) = / f a -V a dfl+ p a -V a dS, (10) 

/• re a (v a ) 

h a {v a ) =3 / Ma / zw„(z) dz dil , (11) 

where p a G |"JJ"nn^ 2 (Sa)1 3 , and f Q G [-^(^a)] 3 is the vector of volume forces. 

Using [12j , we have shown that the original contact problem has an alterna- 
tive weak formulation as the following minimization problem on the set K: 

F(u) = A (u, u) /2 - H(u) - L (u) -> min . (12) 
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Bilinear form A is symmetric, continuous with constant Ma > and coercive 
with constant Ba > 0, and linear form L is continuous. Nonquadratic functional 
H is doubly Gateaux differentiable in Vq: 

H '( U ' y ) = ^2 h 'a( u a, V„), ff''(u,V,w) =^/l£(u a ,V Qt) W Q! ), U,V,W £ V , 
a a 

(13) 

h' a (u a ,v a )=2 fa^ a (e a (u tt )) y^e a ij(u a )e a ij-(v a ) dSl. (14) 

Moreover, we have proved that the following conditions hold: 

(3C>0)(VueF ) {(l-C)A(u,u) >2H{u) }, (15) 

(VueFo)(3i?>0)(VveVb){|fl"(u,v)|<i?||v|| Vb }, (16) 

(3D > 0) (Vu,v,w G V ) { |ff"(u,v,w)| < D \\v\\ Vf> \\w\\ Vo } , (17) 

(3B > 0) (Vu, v G Vo) {A(v,v) - H"(u,v,v) > B ||v|ft o } . (18) 

From these properties, it follows that there exists a unique solution u G K 
of minimization problem (|12[) . and this problem is equivalent to the following 
variational inequality, which is nonlinear in u: 

A(u,v-u)-i?'(u,v-u)-i(v-u) >0, Vveif, ueK. (19) 

To obtain a minimization problem in the whole space Vo, we apply a penalty 
method [1,0,0, El to problem (fT2j) . We use a penalty in the form 



2 

L 2 (S a p) 



{a, P}EQ 

+ Ye ^ \\ u -'- u P'W 2 [L 2 (r a , ,)]s> (2°) 



{a', P'jEQ 1 

where 9 > is a penalty parameter, and y~ = min{0, y}. 

Now, consider the following unconstrained minimization problem in Vo: 

F (u) = A(u,u)/2-JJ(u)-i(u) + J fl (u)-> min . (21) 

u e vb 

The penalty term Jg is nonnegative and Gateaux differentiable in Vo, and its 
differential Jg(u, v) = -\ J2{ a , p} e Q Is afl ( d <*P ~ Uan ~ u 0n)' {v an + vp n ) dS+ 
\ S{ a ' 0'} eq 1 It , , ( u "' _ u P')-(y a > ~ vp) dS satisfy the following properties 
0: 

(VuG7o)(3-R>0)(VvGVb){l4(u,v)| <fl||v|| Vo }, (22) 

(3^>0)(Vu,v,wey ){|4(u + w,v)-^(u,v)|<^||v|| yo ||w|| l/o }, (23) 

(Vu, v G V ) {J' e (u + v, v) - J' e (u, v) > 0} . (24) 

Using these properties and the results in Q, we have shown that problem 
(|2ip has a unique solution u.g G Vo and is equivalent to the following nonlinear 
variational equation in the space Vo: 

^(u,v)=A(u ) v)- J E'(u,v)+4(u,v)-i(v) = 0, VveVo, u G V . (25) 

Using the results of works @, EH, we have proved that ||ue — u||y o g ^ Q ^- 
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4 Iterative methods for nonlinear variational 
equations 

In arbitrary reflexive Banach space Vq consider an abstract nonlinear variational 
equation 

$(u,v) =£(v), VveVb, ueV (26) 

where $ : Vo X Vo — > R is a functional, which is linear in v, but nonlinear in 
u, and L is linear continuous form. Suppose that this variational equation has 
a unique solution u» € Vq. 

For the numerical solution of (|2l)f we use the next iterative method 0, [l5| : 

G(u fe+1 ,v) = G(u fe ,v)- 7 [$ (u fc , v) - L (v)] , Vv e V , k = 0,l,..., (27) 

where G is some given bilinear form in Vq x Vo, 7 £ M. is fixed parameter, and 
u fe e Vo is the fc-th approximation to the exact solution of problem (|2l)|) . 

We have proved the next theorem [^, about the convergence of this 
method. 

Theorem 1. Suppose that the following conditions hold 

(Vu e Vo) > 0) (Vv e Vo) { |$ (u, v)| < \\v\\ Vo }, (28) 
(3D $ >0)(Vu,v,wgV ){|$(u + w ) v)-$(u,v)| < £>*||v|| Vo ||w|| , (29) 
{3B^ > 0) (Vu, v e Vo) { $ (u + v, v) - $ (u, v) > \\vf Vo } , (30) 

bilinear form G is symmetric, continuous with constant Mq > and coercive 
with constant Bq > 0, and 7 € (0; 27*) ,7* = B^Bqj { D\. 

Then ||u fe — u* \\ v — ¥ 0, where {u fe } C Vo is obtained by method (|27[). 

Moreover, the convergence rate in norm \\-\\ G = \/ G (•, •) is linear, and the 
highest convergence rate in this norm reaches as 7 = 7*. 

In addition, we have proposed nonstationary iterative method to solve (l26l) , 
where bilinear form G and parameter 7 arc different in each iteration: 

G fc (u fe+1 ,v) = G fc (u fe ,v) - 7 fe [$(u fe ,v) -L(v)] , Vv e V , k = 0, 1, ... . 

(31) 

A convergence theorem for this method is proved in [15j . 

5 Domain decomposition schemes for contact 
problems 

Now let us apply iterative methods (|2T|) and ([3~T]) to the solution of nonlinear 
penalty variational equation (|25p of multibody contact problem. This penalty 
equation can be written in form (|26l) . where 

$(u,v) = J 4(u,v)-i/'(u,v) + 4(u,v), u,veV . (32) 

We consider such variants of methods (|2"T|) and pip , which lead to the domain 
decomposition. 
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Let us take the bilinear form G in iterative method (j2"7| as follows 0, [l5[ : 



G(u,v) = 4(u,v)+X(u,v), u,vg4 (33) 



1 w 

Q — 1 



^ / u an v an ip a pdS + ^2 / u Q • v Q Q/3 / dS 1 



where^(x) = {l, xe 5^}V{0, xe 5 Qf ,\Si (J }and^(x) = {l, xe T^,}V 
{0, x G r Q /3'\rQ^/ } are characteristic functions of arbitrary subsets S^p C S Q/3 , 
T^, C r Q( g' of possible unilateral and ideal contact areas respectively. 

Introduce a notation u fc+1 = [u fc+1 — u k ]/j+u k <G Vq. Then iterative method 
(|27p with bilinear form (|33p can be written in such way: 

A (u fe+1 , v) + X (u k+1 , v) = L (v) + X (u fe , v) + H' (u k , v) - J' e (u k , v) , (34) 

u fe+1 = 7 u fc+1 + (l-7)u fc , fc = o,l,... . (35) 

Bilinear form X is symmetric, continuous with constant Mx > 0, and non- 
negative fl5j . Due to these properties, and due to the properties of bilinear 
form A, it follows that the conditions of Theorem 1 hold. Therefore, we obtain 
the next proposition: 

Theorem 2. The sequence {u k } of the method {3$ - (ESP converges 
strongly to the solution of penalty variational equation (|25p for 7 S (0; 2B^Bq / D\), 
where Bq — Ba, B$, — B, £>$ = Ma + D + D. The convergence rate in norm 
|| • || G is linear. 

As the common quantities of the subdomains are known from the previous 
iteration, variational equation splits into N separate equations for each 
subdomain VL a , and method (|34[) - (|35l) can be written in the following equivalent 
form: 

a Q (u^,v Q )+ £ / tfu^v^dS+Ysl ^f< +l -v a dS = 



L(v a ) + ^ I ^ a pu k an + (d a p 



u k 



dS- 



Da/3' 



u k a +(u%-u k a )]-v a dS+ti a (u k a> v a ), Vv a e^ (36) 



< +1 =7i« +1 + (l-7)<, a = l,2,...,N, k = 0, 1, ... . (37) 

In each iteration k of method - (|5T)) . we have to solve N linear varia- 
tional equations in parallel, which correspond to some linear elasticity problems 
in subdomains with additional volume forces in Q a , and with Robin boundary 
conditions on contact areas. Therefore, this method refers to parallel Robin- 
Robin type domain decomposition schemes. 

Taking different characteristic functions ip a p and 4> a '/3', we can obtain dif- 
ferent particular cases of penalty domain decomposition method (I56l - ([3T|) . 

Thus, taking 4> a p(x) = 0, (3 £ B a , tf> a pi(x) = 0, [3' € I a , a = 1,2, ...,N, we 
get parallel Neumann-Neumann domain decomposition scheme. 
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Other borderline case is when i/) a p(x) = 1, j3 € B a , {j> a pi(x) = 1, f3' € I a , 
a = 1,2, ...,iV, i.e. = S a p, T 1 ^, = T a p>. 

Moreover, we can choose functions ip a p and c/) a pi differently in each iteration 
k. Then we obtain nonstationary domain decomposition schemes, which are 
equivalent to iterative method (1311) with bilinear forms 



G fc (u,v) = A(u,v) +A fc (u,v), u,veV , fc = 0,l,... , (38) 



a=l 



^ / U an V an ll> af} dS + ^ / u a • v a <f>ap> d S 



/3eB a " °<*l 3 /3>£l 

If we take characteristic functions tp^p and (f)^, as follows d, [3, [H| : 

W X) " *°"> W ~ \ 1, d^(x) - <„(x) - uj|„(x') < ' X - X G ba ? 

^,(x) = l, xeT^, ^€B a! /3'e/ a , a = 1,2, ...,iV, 

then we shall get the method, which can be conventionally named as nonsta- 
tionary parallel Dirichlet-Dirichlet domain decomposition scheme. 

In addition to methods (j2"T)l . (J33J) and (|3"Tj) , (j3"5)) . we have proposed another 
family of DDMs for the solution of (j2"5)l , where the second derivative of functional 
H(u) is used. These domain decomposition methods are obtained from (f3~Tj). if 
we choose bilinear forms G k (u, v) as follows 

G fc (u,v) =A(u,v)-iJ"(u fc ,u,v)+X fe (u,v), u,veV , fc = 0, 1, ... . (39) 

Numerical analysis of presented penalty Robin-Robin DDMs has been made 
for plane unilateral two-body and three-body contact problems of linear elas- 
ticity {u a = 0) using finite element approximations @, E3, Numerical ex- 
periments have confirmed the theoretical results about the convergence of these 
methods. 

Among the positive features of proposed domain decomposition schemes are 
the simplicity of the algorithms and the regularization of original contact prob- 
lem because of the use of penalty method. These domain decomposition schemes 
have only one iteration loop, which deals with domain decomposition, nonlin- 
earity of the stress-strain relationship, and nonlinearity of unilateral contact 
conditions. 
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